#include "Struct_Analysis.h"

#include "Struct_Data.h"

#include "Node_Base.h"

#include "Element_Beam_Twist.h"
#include "Element_Beam3_Twist.h"

#include "Boundary_Base.h"

#include "Load_Node.h"

Struct_Analysis::Struct_Analysis(const Struct_Data& data) :m_data(data)
{

}

bool Struct_Analysis::Analyse()
{
	AssignDOF();
	AssembleKs();
	TreatBoundary();
	AssembleLoad();
	solve();
	Post_Processing();

	return true;
}

void Struct_Analysis::AssignDOF()
{
	for (auto& a : m_data.m_Element)
	{
		auto pElemnet = a.second;
		int nDOFNode = pElemnet->getDOFNode();
		for (auto& b : pElemnet->m_pNode)
		{
			auto pNode = b.lock();
			if (pNode->m_DOF.size() < nDOFNode)
			{
				pNode->m_DOF.resize(nDOFNode);
			}
		}
	}

	for (auto& a : m_data.m_Node)
	{
		auto pNode = a.second;
		for (auto& b : pNode->m_DOF)
		{
			b = -1;
		}
	}

	int iStart = 0;
	for (auto& a : m_data.m_Boundary)
	{
		auto pBoundary = a.second;
		auto pNode = pBoundary->m_pNode.lock();
		int iDirection = pBoundary->m_Direction;
		if (-1 == pNode->m_DOF[iDirection])
		{
			pNode->m_DOF[iDirection] = iStart++;
		}
	}
	m_nFixed = iStart;



	for (auto& a : m_data.m_Node)
	{
		auto pNode = a.second;
		for (auto& b : pNode->m_DOF)
		{
			if (b == -1)
			{
				b = iStart++;
			}
		}
	}
	m_nFree = iStart - m_nFixed;
}

void Struct_Analysis::AssembleKs()
{
	m_K11.resize(m_nFixed, m_nFixed);
	m_K21.resize(m_nFree, m_nFixed);
	m_K22.resize(m_nFree, m_nFree);
	typedef Eigen::Triplet<double> Tri;
	std::list<Tri> L11, L21, L22;

	Eigen::MatrixXd ke;
	std::vector<int> DOFs;
	for (auto& a : m_data.m_Element)
	{
		auto pElement = a.second;
		pElement->get_Ke(ke);
		pElement->getDOFs(DOFs);
		auto nDOF = DOFs.size();
		for (int i = 0; i < nDOF; i++)
		{
			int ii = DOFs[i];
			for (int j = 0; j < nDOF; j++)
			{
				int jj = DOFs[j];
				auto kij = ke(i, j);
				if (ii < m_nFixed && jj < m_nFixed)
				{
					L11.push_back(Tri(ii, jj, kij));
				}
				else if (ii >= m_nFixed && jj < m_nFixed)
				{
					L21.push_back(Tri(ii - m_nFixed, jj, kij));
				}
				else if (ii >= m_nFixed && jj >= m_nFixed)
				{
					L22.push_back(Tri(ii - m_nFixed, jj - m_nFixed, kij));
				}
			}
		}
	}
	m_K11.setFromTriplets(L11.begin(), L11.end());
	m_K21.setFromTriplets(L21.begin(), L21.end());
	m_K22.setFromTriplets(L22.begin(), L22.end());
}

void Struct_Analysis::TreatBoundary()
{
	x1.resize(m_nFixed);
	for (auto& a : m_data.m_Boundary)
	{
		auto pBoundary = a.second;
		auto pNode = pBoundary->m_pNode.lock();
		auto DOF = pNode->m_DOF[pBoundary->m_Direction];
		x1[DOF] = pBoundary->m_Value;
	}
}

void Struct_Analysis::AssembleLoad()
{
	F1.resize(m_nFixed);
	F2.resize(m_nFree);
	F1.setZero();
	F2.setZero();
	for (auto& a : m_data.m_Load)
	{
		auto pLoad = a.second;
		pLoad->Assemble(F1, F2);
	}
}

void Struct_Analysis::solve()
{
	Eigen::SimplicialLDLT<SpMat> solver;
	solver.analyzePattern(m_K22);
	solver.factorize(m_K22);
	x2 = solver.solve(F2 - m_K21 * x1);
	std::cout << m_K22 << std::endl;
	std::cout << x1 << std::endl << x2 << std::endl;
}

void Struct_Analysis::Post_Processing()
{
	for (auto& a : m_data.m_Element)
	{
		auto& pEntity = a.second;
		pEntity->getparameterlist(x1, x2);
	}
}
